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Abstract: This work concerns causal acoustical wave equations which imply frequency power- 
law attenuation. A connection between the five-parameter fractional Zcncr wave equation, which 
is derived from a fractional stress-strain relation plus conservations of mass and momentum, 
and the physically well established multiple relaxation framework is developed. It is shown that 
for a certain continuous distribution of relaxation mechanisms, the two descriptions are equal. 

Keywords: Wave equations, attenuation, multiple relaxation, fractional modeling, power law 
descriptions, frequency dispersion. 
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1. INTRODUCTION 

The multiple relaxation mechanism framework of Nach- 
man et al. (1990) is widely considered as adequate for 
acoustic wave modeling in lossy complex media like those 
encountered in medical ultrasound. It relics on thermo- 
dynamics and first principles of acoustical physics. The 
corresponding wave equation for N relaxation mechanisms 
is a causal partial differential equation with its highest 
time derivative order N + 2. We denote this the Nachman- 
Smith-Waag (NSW) model. 

Attenuation in complex media often follows a power law: 
afe(u;) cx w'', with rj e [0.2] (Szabo and Wu (2000)). 
The range where experiments indicate this may cover 
many frequency decades. In order to make the NSW 
model attenuation adequately follow lu'^, either the valid 
frequency band must be narrow, or the number of assumed 
mechanisms N must be large thus inferring a partial 
differential equation of very high order. 

Another way to derive a lossy wave equation is to com- 
bine the principles of mass and momentum conservation 
with some stress strain relation. This constitutive relation 
may include fractional time-derivatives, exemplified by the 
fractional Zener model by Holm and Nasholm (2011). The 
resulting wave equation is causal and the corresponding at- 
tenuation follows power laws within wide frequency bands. 

The purpose of the present work is to demonstrate the 
link between the NSW and the fractional Zener models via 
a continuum of relaxation mechanisms. Relevant parts of 
Nasholm and Holm (2011) are reviewed and rcformiilated. 
In addition a more general deduction is provided where 
the fractional Zener model parameters a and ^ are not 
necessarily equal. We aim to encourage the acoustical 
community to more frequently adopt fractional calculus 
descriptions for wave modeling in complex media. 

* This research was partly supported by the "High Resolution 
Imaging and Beamforming" project of the Norwegian Research 
Council. 



2. THEORY 

2.1 Conservation laws and Generalized Compressibility 

The linearized conservation of mass corresponds to the 
strain being defined by 

J- 

e{t) = Vu{x,t),< > e(w) = -iku{k,w), (1) 

where u is the displacement and the symbol J" denotes 

transformation into the spatio-temporal frequency domain 
where lo is the angular frequency and k the wavcmimbcr. 

The linearized conservation of momentum is expressed as 



Va(t) = po 



d^u{x,t) J- 



dt^ 



> —ika{ijj) = po{iu})'^u{k,Lj), 



(2) 



where po is the steady-state mass density and a denotes 
the stress, which in this context corresponds to the nega- 
tive of the pressure. 

The frequency-domain generalized compressibility is de- 
fined as the ratio between strain and stress: k{oj) = 
€{uj)/a{u)), therefore being related to the constitutive 
stress-strain relation. Combining this definition with the 
conservation laws (1) and (2) gives 



k'^{u}) — uj'^pok{uj) 



^ V^u{x,t) - 



dt2 



K{t) * u{x,t) 



= 0. 



(3) 
(4) 



Under circumstances where the linearized conservations of 
mass and momentum are valid, the wave equation is thus 
completely determined by the generalized compressibility. 

The generalized compressibility k{uj) as given above is 
sometimes (e.g. in viscoelasticity) called complex com- 
pliance J*{ui) = 1/G*{(jj), where G*(w) is the complex 
modulus. 



2.2 A Continuum of NSW Relaxation Processes 

The NSW model of multiple discrete relaxation processes 
results in the generalized compressibility 

JV 

«(a;) = «o-iu;^^^^, (5) 



where the mechanisms u = 1 . . . N, have the relaxation 
times Ti , . . . , Tjv and the compressibility contributions 
Ki,. . . ,kn (Nachman et al. (1990)). 

Following Nasholm and Holm (2011), a representation of 
(5) when considering a continuum of relaxation mecha- 
nisms distributed in the frequency band e [01,^2] 
with the compressibility contributions described by the 
distribution K^(f2) becomes 

Kn{uj) = Ko — iuj I '^^^) 

Letting the limits of the integral go between = 

and ^2 = 00, and instead incorporating any possible 
relaxation distribution bandwidth limitation into Avi/(f2), 
the integral above is a Stieltjes transform. Applying the 
Laplace transform relation 

the generalized compressibility (6) becomes 

POO roo 

kn(w) = Ko-iuj / K^{n) / e-"*e-'"*di dfl 
Jo Jo 

= Ko- iLoJ't{H{t)Cn {«;,(r!)Xt)}(w). (8) 
2.3 The Fractional Zener Wave Equation 

The five-parameter fractional Zener stress-strain constitu- 
tive relation is experimentally shown to be valid for a wide 
range of complex media, see the references in Holm and 
Nasholm (2011). As given by Bagley and Torvik (1983), it 
may be expressed as 



a{t) 



dtf^ 



e{t) 



.9"e(i) 



dt<^ 



(9) 



From this relation, the frequency-domain fractional Zener 
compressibility is obtained through the ratio e{w)/cr{(jo): 

Kz(w) - - 



Kq 



1 + {Taioj)"' 

{iuy 



(10) 

Due to thermodynamic constraints, /3 is restricted to be 
smaller than or equal to a (Glockle and Nonnenmacher 

(1991)). 

Insertion of the generalized compressibility (10) into the 
dispersion relation (4), generates the time-domain frac- 
tional Zener wave equation (Holm and Nasholm (2011)) 

1-2. 



2 



d dtP+^ 



= 0. 



(11) 



2.4 Connecting the NSW and the Fractional Zener Models 

Provided that the conservations of mass (1) and momen- 
tum (2) are valid, and provided that the NSW generalized 



compressibility k-^(oj) of (8) is equal to the fractional Zener 
generalized compressibility kz{uj) of (10), the dispersion 
relations from (3) are also equal. Because the dispersion 
relation is a spatio-temporal Fourier representation of the 
wave equation, kn(w) = kz{u)) thus implies that the NSW 
wave equation becomes equal to the fractional Zener wave 
equation (11). Direct comparison of kn(w) in (8) to k-z{oj) 
in (10), tells that they are equal in case the following is 
true: 



J^t{H{t)£n{KAmit)}{^ 



Ko 



(rf/r,")(za;) 



a-{a-fi + l) 



-. (12) 



First we choose to study the case a = (3, which was 
also treated in Nasholm and Holm (2011). Inverse Fourier 
transformation of both sides of (12), then gives 

H{t)Cn{^i,Am(t) 



{iuj) 



a-l 



= Ko{l - T'^/T^)H{t)E^^l {--{t/T„r) , (13) 

where Ea,b{-) is the Mittag-LefFler function (see Appendix 
A), and H{t) is the Heaviside step function. The Fourier 
transform relation used in the last step above is given in 
(A. 2). Moreover, the inverse Laplace transform relation of 
(A. 3), Eq. (13) hence gives 
K,{n) = Ko(l - r^lT:)U,i {n,T-") 

_ 1 Ko(t« - r,")17"-i sin(a7r) 



TT (r<^0)2« + 2(t^Q)« cos(a7r) -|- 1 



'«l/Ml(^) 

(14) 

where /„_i($l, a) was inserted from (A. 4). 

For the more general case /3 < a, inverse Fourier transform 
on both sides of (12) instead gives 
H{t)Cn{K,m]{t) = 

(iu;)"-i 



i Ta + (iw)" J 

Proceeding in a similar manner as for the a = ji case then 
gives the distribution 

«.(n) = KoU,i (f^,T-") - Ko(rf/r«)/„,a_^+i (n,r-«) 

_ kqt^ sin(a7r) 

~ TT (t„1^)2" + 2(r<^f])" cos(a7r) + 1 

Atorf il^-^ sin(^7r) - r^fi^+'^-i sin((a - /3)7r) 



TT 

= '^tML(^)- 



(r^0)2« + 2{t^Q.Y cos(a7r) + 1 



(16) 



We have thus shown that the fractional Zener wave equa- 
tion (11) may be obtained within the Nachman-Smith- 
Waag framework of multiple relaxation (Nachman et al. 
(1990)), when assuming a continuum of relaxation mech- 
anisms with the compressibility contribution as given by 
the distribution k^]vil(^) (l^)- 

In viscoelasticity, a relaxation time spectrum (below de- 
noted H{t)) related to K^iP) is commonly studied (see e.g. 
Glockle and Nonnenmacher (1991) and references therein). 
It is related to the complex modulus through 
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/oo 
fl-(r)e-*/Mlnr. (17) 
-co 

It may be shown that for the 5-parameter fractional Zener 
model, when setting fl = t~^, the r-depcndcncy of H{t) 
differs by a factor r to Ki,ml{^) of (16). Figs. 5 and 
6 of Glockle and Nonnenmacher (1991) illustrate that 
a = /3 gives symmetric H{t), while a ^ f3 breaks the 
symmetry, most significantly far from the peak region. 
Such relaxation spectra are experimentally observed for 
complex media, e.g. natural rubber. 

3. ATTENUATION AND PHASE VELOCITY 
EXAMPLES 

The conventional decomposition of the frequency-dependent 
wavenumber into its real and imaginary parts, gives the 
phase velocity Cp{oj) = co/^{k} and attenuation ak{u}) = 

In general, the attenuation and the phase velocity are thus 
given from the dispersion relation (3) as 

1/2 } r—i (18) 

For the fractional Zener wave equation, this results in three 
distinct regions with attenuation power-laws (Nasholm 
and Holm (2011)): ak oc w^"'"" in a low-frequency regime, 
ak oc in an intermediate frequency regime, and 

ak oc in a high-frequency regime. 

In the following, the fractional Zener phase velocities and 
attenuations are further investigated numerically for the 
a = j3 case in a similar manner as in Nasholm and Holm 
(2011). This is done explicitly by insertion of k{oj) = kn(w) 
into (18). The results from such calculations are compared 
to what is found by insertion of the distribution K,yMh{^) 
of (14) into the NSW generalized compressibility integral 
formula (6). This generalized compressibility is finally 
applied to (18), from which afe(w) and Cp(w) are found. 

We use the latter calculation method to explore the 
effect of letting the continuum of relaxation mechanisms 
populate only a bounded frequency interval, rather than 

the entire G [0, oo] region. 

Figure 1 compares attenuation curves and shows the dis- 
tributions K^{il), while Fig. 2 displays the correspond- 
ing frequency-dependent phase velocity. Note the high- 
frequency asymptote of K^{fl), which has the fractal prop- 
erty of being proportional to The integral over 
n G [r2i,r22] in the calculation of kn(w) from (6) is eval- 
uated numerically using the recursive adaptive Simpson 
quadrature method. 

4. CONCLUSION 

This work shows analytically that the lossy fractional 
Zener wave equation (9) (Holm and Nasholm (2011)) 
may be attained within the NSW multiple relaxation loss 
framework (Nachman et al. (1990)), given that a contin- 
uum of relaxation processes are weighted appropriately 
following KyML(^) described in (16). The result may be 
seen as a generalization of what was presented in Nasholm 



and Holm (2011) as the developments here also cover the 
case of the fractional Zener constitutive relation not having 
equal derivative orders a and /?. 

The advantages of the fractional Zener model is that it 
fits measurements well and that it is characterized by a 
small number of parameters, while the NSW model is more 
intuitive as it does not comprise fractional derivatives. It 
is also better rooted in fundamental physics. 

REFERENCES 

Bagley, R.L. and Torvik, P.J. (1983). Fractional calculus 
— A different approach to the analysis of viscoelastically 
damped structures. AIAA J., 21(5), 741-748. 

Djrbashian, M.M. (1966). Integral transforms and repre- 
sentations of functions in the complex domain, chapter 
3-4. Nauka, Moscow, USSR. In Russian. 

Djrbashian, M.M. (1993). Harmonic analysis and bound- 
ary value problems in the complex domain, chapter 1. 
Birkhauser, Basel, Switzerland. 

Glockle, W.G. and Nonnenmacher, T.F. (1991). Fractional 
integral operators and Fox functions in the theory of 
viscoelasticity. Macromolecules, 24(24), 6426-6434. 

Haubold, H.J., Mathai, A.M., and Saxena, R.K. (2011). 
Mittag-LefFler functions and their applications. Journal 
of Applied Mathematics, 2011, 1 51. 

Holm, S. and Nasholm, S.P. (2011). A causal and fractional 
all-frequency wave equation for lossy media. J. Acoust. 
Soc. Am., 130(4), 2195-2202. 

Mittag-Leffer, M.G. (1903). Sur la nouvelle fonction Ea{x) 
(On the new function Ea{x)). C. R. Acad. Sci. Paris, 
137, 554-558. 

Nachman, A.I., Smith III, J.F., and Waag, R.C. (1990). 
An equation for acoustic propagation in inhomogeneous 
media with relaxation losses. J. Acoust. Soc. Am., 88, 
1584-1595. 

Nasholm, S.P. and Holm, S. (2011). Linking multiple 
relaxation, power-law attenuation, and fractional wave 
equations. J. Acoust. Soc. Am., 130(5), 3038-3045. 

Podlubny, I. (1999). Fractional differential equations, 
chapter 1-2. Academic Press, New York. 

Szabo, T.L. and Wu, J. (2000). A model for longitudinal 
and shear wave propagation in viscoelastic media. J. 
Acoust. Soc. Am., 107, 2437-2446. 

Wiman, A. (1905). Uber den Fundamentalsatz in der 
Theorie der Funktionen Ea{x) (About the fundamental 
theorem in the theory of the function Ea{x)). Acta 
Mathematica, 29, 191-201. 

Appendix A. THE MITTAG-LEFFLER FUNCTION 
A.l Definition and Fourier Transform Relation 

The one-parameter Mittag-Leffler function was introduced 

in Mittag-Leffer (1903). A two-parameter analogy was 
presented in Wiman (1905), which may be written as 

n— 

where F is the Euler Gamma function and the parameters 
are commonly restricted to {a, 6} e C, K{a, 6} > 0, and 
f e C. See Haubold et al. (2011) for a comprehensive 
review of Mittag-LefHer function properties. 



Pre-print paper. Presented at the 5th IFAC Symposium on Fractional Differentiation cind Its 
Applications (FDA 2012), Hohai University, Nanjing, China, 14-17 May 2012 




10"' 



T ■ ^2e[10"^ 10^] 



la „, , T ■^26[10"^ 10^] 

k,ML a ^ ^ 

,a „, , T ■£26[10"^ 10^] 

k,ML a ^ ^ 



10" 



10" 



10" 



10" 







T =1000t , a=0.3 

a E 




vML a 


£2e[10"^, 10^]^^^^,^ 




vML a 


liellO-^, 10^] 




;- vML a 


i2e[10"'', 10^] 




10-= 


10" 

o 


10= 




10" 



'a, T £26 [10"^, 10^] 

a, T ■£2E[10-^ 10^] 

k,ML a ^ 

,a, T ■£2e[10-^ 10^] 

k,ML a ^ 



y 10" 



10-= 


10° 

COT 

a 


10= 






T =1000x , a=0.8 

a e 






X ■ ile[10"^ 10^1 




- ^— ^ML' 


X ■ ^le[10"^ 10^1 




- -;- \mU 


X ■ ile[10"^ 10^1 

a ; 




10-= 


10° 

a 


10= 



Fig. 1. Top panes: frequency-dependent attenuation for To- = lOOOre with the fractional derivative orders a = 0.3 
(top-left pane), and a — 0.8 (top-right pane). The attenuation curves display both explicit calculations from the 
fractional Zener model, and calculations using the distribution k^ml{^) in the NSW model. Different choices of 
integration limits for fii and in (6) are exploited, as displayed in the legends. The horizontal axis represents 
normalized frequency. For visualization convenience, each attenuation curve is normalized to = 1 at ujTo- = 1. 
Bottom panes: the corresponding normalized effective compressibilities k,/ml(^) of the continuum of relaxation 
processes as a function of normalized relaxation frequency ft ■ for a = 0.3 (bottom-left pane), and a = 0.8 
(bottom- right pane). The £7 integration limits are given in the legends. 
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Fig. 2. Frequency-dependent phase velocity for Ta = IOOOte with the fractional derivative orders a = 0.3 (left pane), 
and a = 0.8 (right pane). The curves display Cp{(jj) as predicted by the fractional Zener model, as well as predicted 
by using the distribution K^ML(f^) in the NSW model. Different choices of integration limits for Vli and VI2 in (6) 
are exploited, as displayed in the legends. The horizontal axis represents normalized frequency. 

A useful Fourier transform pair involving the Mittag- ,_ (10;)"^^^ 

LefBer function is (Podlubny (1999)) ^ {H{t) t Ea,b{ ~ At") } (co) =^^:j—^- (A.2) 
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A. 2 Laplace Transform Integral Representation 

The function t''^'^ Ea.b{-At°-) may for < a < 1 be written 
on an integral form (Djrbashian (1966, 1993)): 



poo 

t''-^E„4-Af')= e-^Va,6(f^,^)dO, (A.3) 

^0 



where 



f (O - -4sin[(b - g)^] + »° sin(b7r) 

Ja,bK^L,A)- ^ 02a + 2A0«C0s(a7r) + A2 ■ ^"^"^^ 
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